---
title: "R Notebook"
output: html_notebook
---

This is an [R Markdown](http://rmarkdown.rstudio.com) Notebook. When you execute code within the notebook, the results appear beneath the code. 

Try executing this chunk by clicking the *Run* button within the chunk or by placing your cursor inside it and pressing *Cmd+Shift+Enter*. 

```{r}
library(survival)
library(tidyverse)
library(survminer)
library(texreg)
library(stargazer)
library(MatchIt)
library(cobalt)
```


```{r dataset}
cf = read.csv("cf.csv")
legitimacy = read.csv("legitimacy.csv")
pa = read.csv("pa.csv")
```

```{r ceasefires}

cf$catcoup = factor(cf$catcoup, levels = c("no coup", "failed", "success"))
cf$catcoupm = factor(cf$catcoupm, levels = c("no coup", "failed", "success"))

cf1 = coxph(Surv(istart_onset, ifinish_onset, all_cf) ~ coup + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num + terrcont + strata(gcf_collap) + frailty(gwno_a,distribution="gamma",method="em"),
                     method=c("efron"), data = cf)


cf2 = coxph(Surv(istart_onset, ifinish_onset, all_cf) ~ coup + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terrcont + multirebel + territorial + strata(gcf_collap) + frailty(gwno_a,distribution="gamma",method="em"),
                     method=c("efron"), data = cf)


cf3 = coxph(Surv(istart_onset, ifinish_onset, all_cf) ~ catcoup + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num + terrcont + strata(gcf_collap) + frailty(gwno_a,distribution="gamma",method="em"),
                     method=c("efron"), data = cf)


cf4 = coxph(Surv(istart_onset, ifinish_onset, all_cf) ~ catcoup + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terrcont + multirebel + territorial + strata(gcf_collap) + frailty(gwno_a,distribution="gamma",method="em"),
                     method=c("efron"), data = cf)

cf1h = exp(coef(cf1))
cf2h = exp(coef(cf2))
cf3h =  exp(coef(cf3))
cf4h = exp(coef(cf4))

texreg(list(cf1, cf2, cf3, cf4), include.aic = FALSE, include.rsquared = FALSE, include.zph = FALSE, include.missings = FALSE, override.coef = list(cf1h, cf2h, cf3h, cf4h), custom.coef.names = c("Coup", "GDP (logged)", "Level of Democracy", "Intensity Level", "Conflict Duration (logged)", "Relative Strengths", "Territorial Control", "Multirebel Conflict", "Territorial Conflict", "Failed Coup", "Successful Coup"))


### Military Coup ONLY and Ceasefires

cf5 = coxph(Surv(istart_onset, ifinish_onset, all_cf) ~ military_coup + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num + terrcont + strata(gcf_collap) + frailty(gwno_a,distribution="gamma",method="em"),
                     method=c("efron"), data = cf)


cf6 = coxph(Surv(istart_onset, ifinish_onset, all_cf) ~ military_coup + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terrcont + multirebel + territorial + strata(gcf_collap) + frailty(gwno_a,distribution="gamma",method="em"),
                     method=c("efron"), data = cf)


cf7 = coxph(Surv(istart_onset, ifinish_onset, all_cf) ~ catcoupm + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num + terrcont + strata(gcf_collap) + frailty(gwno_a,distribution="gamma",method="em"),
                     method=c("efron"), data = cf)


cf8 = coxph(Surv(istart_onset, ifinish_onset, all_cf) ~ catcoupm + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terrcont + multirebel + territorial + strata(gcf_collap) + frailty(gwno_a,distribution="gamma",method="em"),
                     method=c("efron"), data = cf)


cf5h = exp(coef(cf5))
cf6h = exp(coef(cf6))
cf7h =  exp(coef(cf7))
cf8h = exp(coef(cf8))



texreg(list(cf5, cf6, cf7, cf8), include.aic = FALSE, include.rsquared = FALSE, include.zph = FALSE, include.missings = FALSE, override.coef = list(cf5h, cf6h, cf7h, cf8h), custom.coef.names = c("Military Coup", "GDP (logged)", "Level of Democracy", "Intensity Level", "Conflict Duration (logged)", "Relative Strengths", "Territorial Control", "Multirebel Conflict", "Territorial Conflict", "Failed Military Coup", "Successful Military Coup"))


sum_cf = cf %>% mutate(log_gdp = log(gdp), log_time = log(time)) %>%  select(all_cf, coup, military_coup, log_gdp, v2x_polyarchy, intensity_level, log_time, rebstrength_num, terrcont, multirebel, territorial)

stargazer(as.data.frame(sum_cf), covariate.labels = c("Ceasefire", "Coup", "Military Coup", "GDP (logged)", "Level of Democracy", "Intensity Level", "Conflict Duration (logged)", "Relative Strengths", "Territorial Control", "Multirebel Conflict", "Territorial Conflict"))
```

```{r}
mod_cf = coxph(Surv(istart_onset, ifinish_onset, all_cf) ~ coup + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terrcont + multirebel + territorial, data = cf)

# Create new data for prediction
new_cf = data.frame(
  coup = c(0, 1),  
  gdp = mean(cf$gdp, na.rm = TRUE),
  v2x_polyarchy = mean(cf$v2x_polyarchy, na.rm = TRUE),
  intensity_level = mean(cf$intensity_level, na.rm = TRUE),
  time = mean(cf$time, na.rm = TRUE),
  rebstrength_num = mean(cf$rebstrength_num, na.rm = TRUE),
  terrcont = mean(cf$terrcont, na.rm = TRUE),
  multirebel = mean(cf$multirebel, na.rm = TRUE),
  territorial = mean(cf$territorial, na.rm = TRUE)
)


surv_cf = survfit(mod_cf, newdata = new_cf)

ggsurvplot(surv_cf, data = cf, 
           conf.int = FALSE,   
           legend.title = "", 
           legend.labs = c("No Coup", "Coup"),
           title = "",
           xlab = "Time (Months)",
           ylab = "Survival Probability",
           ggtheme = theme_minimal(), 
           palette = c("lightgrey", "darkgray"),  
           linetype = c("dotted", "solid"))



```






```{r Matching}
m.cf1 = matchit(coup ~ log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terrcont + multirebel + territorial, data = cf, method = "full", caliper = 0.01)
matched_cf1 = match.data(m.cf1)

cf_matching1 = coxph(Surv(istart_onset, ifinish_onset, all_cf) ~ coup + log(gdp) + v2x_polyarchy +  
      intensity_level + log(time) + rebstrength_num + terrcont + 
      multirebel + territorial, 
      weights = weights, 
      cluster = gwno_a, 
      data = matched_cf1)


cf_matching2 = coxph(Surv(istart_onset, ifinish_onset, all_cf) ~ coup + log(gdp) + v2x_polyarchy +  
      intensity_level + log(time) + rebstrength_num + terrcont + 
      multirebel + territorial, 
      weights = weights, 
      cluster = subclass, 
      data = matched_cf1)

m.cf2 = matchit(coup ~ log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terrcont + multirebel + territorial, data = cf, method = "full", caliper = 0.5)
matched_cf2 = match.data(m.cf2)



cf_matching3 = coxph(Surv(istart_onset, ifinish_onset, all_cf) ~ coup + log(gdp) + v2x_polyarchy +  
      intensity_level + log(time) + rebstrength_num + terrcont + 
      multirebel + territorial, 
      weights = weights, 
      cluster = gwno_a, 
      data = matched_cf2)


cf_matching4 = coxph(Surv(istart_onset, ifinish_onset, all_cf) ~ coup + log(gdp) + v2x_polyarchy +  
      intensity_level + log(time) + rebstrength_num + terrcont + 
      multirebel + territorial, 
      weights = weights, 
      cluster = subclass, 
      data = matched_cf2)


new.names = c(
  "coup" = "Coup", 
  "log(gdp)" = "GDP (logged)", 
  "v2x_polyarchy" = "Level of Democracy",
  "intensity_level_2" = "Intensity Level",
  "log(time)" = "Conflict Duration (logged)", 
  "rebstrength_num" = "Relative Strengths",
  "terrcont" = "Territorial Control",
  "multirebel" = "Multirebel Conflict",
  "territorial" = "Territorial Conflict",
  "distance" = "Distance"
)

love.plot(m.cf1, var.names = new.names)
love.plot(m.cf2, var.names = new.names)

cfmatching1h = exp(coef(cf_matching1))
cfmatching2h = exp(coef(cf_matching2))
cfmatching3h =  exp(coef(cf_matching3))
cfmatching4h = exp(coef(cf_matching4))



texreg(list(cf_matching1, cf_matching2, cf_matching3, cf_matching4), include.aic = FALSE, include.rsquared = FALSE, include.zph = FALSE, include.missings = FALSE, override.coef = list(cfmatching1h, cfmatching2h, cfmatching3h, cfmatching4h), custom.coef.names = c("Coup", "GDP (logged)", "Level of Democracy", "Intensity Level", "Conflict Duration (logged)", "Relative Strengths", "Territorial Control", "Multirebel Conflict", "Territorial Conflict"))



print(xtable(as.data.frame(bal.tab(m.cf1, var.names = new.names)$Balance), digits = 3), sanitize.rownames.function = function(x) new.names[x])

print(xtable(as.data.frame(bal.tab(m.cf2, var.names = new.names)$Balance), digits = 3), sanitize.rownames.function = function(x) new.names[x])



```



```{r legitimacy-seeking rebels}
cf9 = coxph(Surv(istart_onset, ifinish_onset, cf) ~ coup*healthedu + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) + strata(cf_collap) + frailty(dyad_id, distribution="gamma",method="em"), method=c("efron"), data = legitimacy)


cf9h =  exp(coef(cf9))

cf10 = coxph(Surv(istart_onset, ifinish_onset, cf) ~ coup*healthedu + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terrcont + territorial + strata(cf_collap) + frailty(dyad_id, distribution="gamma",method="em"), method=c("efron"), data = legitimacy)


cf10h =  exp(coef(cf10))


texreg(list(cf9, cf10), include.aic = FALSE, include.rsquared = FALSE, include.zph = FALSE, include.missings = FALSE, override.coef = list(cf9h, cf10h), custom.coef.names = c("Coup", "Service", "GDP (logged)", "Level of Democracy", "Intensity Level", "Conflict Duration (logged)",  "Coup X Service", "Relative Strengths", "Territorial Control", "Territorial Conflict"))


cf11 = coxph(Surv(istart_onset, ifinish_onset, cf) ~ military_coup*healthedu + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) + strata(cf_collap) + frailty(dyad_id, distribution="gamma",method="em"), method=c("efron"), data = legitimacy)
  
  
  
cf12 = coxph(Surv(istart_onset, ifinish_onset, cf) ~ military_coup*healthedu + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terrcont + territorial + strata(cf_collap) + frailty(dyad_id, distribution="gamma",method="em"), method=c("efron"), data = legitimacy)


cf11h =  exp(coef(cf11))
cf12h =  exp(coef(cf12))


texreg(list(cf11, cf12), include.aic = FALSE, include.rsquared = FALSE, include.zph = FALSE, include.missings = FALSE, override.coef = list(cf11h, cf12h), custom.coef.names = c("Military Coup", "Service", "GDP (logged)", "Level of Democracy", "Intensity Level", "Conflict Duration (logged)",  "Military Coup X Service", "Relative Strengths", "Territorial Control", "Territorial Conflict"))


sum_legitimacy = legitimacy %>% mutate(log_gdp = log(gdp), log_time = log(time)) %>%  select(cf, coup, military_coup, healthedu, log_gdp, v2x_polyarchy, intensity_level, log_time, rebstrength_num, terrcont, territorial)

stargazer(as.data.frame(sum_legitimacy), covariate.labels = c("Ceasefire", "Coup", "Military Coup", "Service", "GDP (logged)", "Level of Democracy", "Intensity Level", "Conflict Duration (logged)", "Relative Strengths", "Territorial Control", "Territorial Conflict"))

```


```{r legitimacy rebel}
legitimacy2 = legitimacy %>%
  mutate(group = case_when(
    coup == 1 & healthedu == 1 ~ "Coup & Service",
    coup == 1 & healthedu == 0 ~ "Coup & No Service",
    coup == 0 & healthedu == 1 ~ "No Coup & Service",
    TRUE ~ "No Coup & No Service"
  ))

legitimacy2$group = factor(legitimacy2$group, 
                            levels = c("No Coup & No Service", "No Coup & Service", 
                                       "Coup & No Service", "Coup & Service"))

mod_legitimacy = coxph(Surv(istart_onset, ifinish_onset, cf) ~ 
                          group + log(gdp) + v2x_polyarchy + 
                          intensity_level + log(time) + rebstrength_num + 
                          terrcont + territorial, 
                        data = legitimacy2)


new_legitimacy = data.frame(
  group = unique(legitimacy2$group),  
  gdp = mean(legitimacy2$gdp, na.rm = TRUE),
  v2x_polyarchy = mean(legitimacy2$v2x_polyarchy, na.rm = TRUE),
  intensity_level = mean(legitimacy2$intensity_level, na.rm = TRUE),
  time = mean(legitimacy2$time, na.rm = TRUE),
  rebstrength_num = mean(legitimacy2$rebstrength_num, na.rm = TRUE),
  terrcont = mean(legitimacy2$terrcont, na.rm = TRUE),
  territorial = mean(legitimacy2$territorial, na.rm = TRUE)
)

surv_legitimacy = survfit(mod_legitimacy, newdata = new_legitimacy)


ggsurvplot(surv_legitimacy, data = new_legitimacy, 
           conf.int = FALSE,   
           legend.title = "", legend.labs = c("No Coup, No Service", 
                           "No Coup, Service", 
                           "Coup, No Service", 
                           "Coup, Service"),  # Custom labels
           ggtheme = theme_minimal(),
           palette = c("grey", "darkgray", "CadetBlue","MidnightBlue"),  # Muted colors
           linetype = "strata",
           size = 0.8,
           xlim = c(0, 100),
           xlab = "Time (Months)") ## It implies that health and education reforms do not significantly alter the probability of a ceasefire in the absence of a coup.

```






```{r} 

m.legitimacy1 = matchit(coup ~ log(gdp) + healthedu + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terrcont +  territorial, data = legitimacy, method = "full", caliper = 0.01, std.caliper = FALSE)
matched_legitimacy1 = match.data(m.legitimacy1)


legitimacy_matching1 = coxph(Surv(istart_onset, ifinish_onset, cf) ~ coup*healthedu + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terrcont + territorial, weights = weights, 
      cluster = dyad_id, data = matched_legitimacy1)



legitimacy_matching2 = coxph(Surv(istart_onset, ifinish_onset, cf) ~ coup*healthedu + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terrcont + territorial, weights = weights, 
      cluster = subclass, data = matched_legitimacy1)



m.legitimacy2 = matchit(coup ~ log(gdp) + healthedu + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terrcont +  territorial, data = legitimacy, method = "full", caliper = 0.02, mahvars = "gdp")
matched_legitimacy2 = match.data(m.legitimacy2)

legitimacy_matching3 = coxph(Surv(istart_onset, ifinish_onset, cf) ~ coup*healthedu + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terrcont + territorial, weights = weights, 
      cluster = dyad_id, data = matched_legitimacy2)


legitimacy_matching4 = coxph(Surv(istart_onset, ifinish_onset, cf) ~ coup*healthedu + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terrcont + territorial, weights = weights, 
      cluster = subclass, data = matched_legitimacy2)


new.names.legit = c(
  "coup" = "Coup", 
  "healthedu" = "Service",
  "log(gdp)" = "GDP (logged)", 
  "v2x_polyarchy" = "Level of Democracy",
  "intensity_level_2" = "Intensity Level",
  "log(time)" = "Conflict Duration (logged)", 
  "rebstrength_num" = "Relative Strengths",
  "terrcont" = "Territorial Control",
  "territorial" = "Territorial Conflict",
  "distance" = "Distance",
  "gdp" = "GDP"
)

legitimacy1h =  exp(coef(legitimacy_matching1))
legitimacy2h =  exp(coef(legitimacy_matching2))

legitimacy3h =  exp(coef(legitimacy_matching3))
legitimacy4h =  exp(coef(legitimacy_matching4))

texreg(list(legitimacy_matching1, legitimacy_matching2, legitimacy_matching3, legitimacy_matching4), include.aic = FALSE, include.rsquared = FALSE, include.zph = FALSE, include.missings = FALSE, override.coef = list(legitimacy1h, legitimacy2h, legitimacy3h, legitimacy4h), custom.coef.names = c("Coup", "Service", "GDP (logged)", "Level of Democracy", "Intensity Level", "Conflict Duration (logged)", "Relative Strengths", "Territorial Control", "Territorial Conflict", "Coup X Service"))

love.plot(m.legitimacy1, var.names = new.names.legit)
love.plot(m.legitimacy2, var.names = new.names.legit)        



print(xtable(as.data.frame(bal.tab(m.legitimacy1, var.names = new.names.legit)$Balance), digits = 3), sanitize.rownames.function = function(x) new.names.legit[x])

print(xtable(as.data.frame(bal.tab(m.legitimacy2, var.names = new.names.legit)$Balance), digits = 3), sanitize.rownames.function = function(x) new.names.legit[x])



```




```{r pa}
#### All Coups ####
pa$catcoup = factor(pa$catcoup, levels = c("no coup", "failed", "success")) 
pa1 = coxph(Surv(istart_pa, ifinish_pa, pa) ~ coup + all_cf + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terr + strata(enum_pa) + frailty(gwno_a,distribution="gamma",method="em"), method=c("efron"), data = pa)

pa2 = coxph(Surv(istart_pa, ifinish_pa, pa) ~ coup + all_cf + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terr + multirebel + territorial  + strata(enum_pa) + frailty(gwno_a,distribution="gamma",method="em"), method=c("efron"), data = pa)


pa3 = coxph(Surv(istart_pa, ifinish_pa, pa) ~ catcoup + all_cf + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terr + strata(enum_pa) + frailty(gwno_a,distribution="gamma",method="em"), method=c("efron"), data = pa)

pa4 = coxph(Surv(istart_pa, ifinish_pa, pa) ~ catcoup + all_cf + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terr + multirebel + territorial  + strata(enum_pa) + frailty(gwno_a,distribution="gamma",method="em"), method=c("efron"), data = pa)

# hazard ratios
pa1h = exp(coef(pa1))
pa2h = exp(coef(pa2))
pa3h =  exp(coef(pa3))
pa4h = exp(coef(pa4))




texreg(list(pa1, pa2, pa3, pa4), override.coef = list(pa1h, pa2h, pa3h, pa4h), include.aic = FALSE, include.rsquared = FALSE, include.zph = FALSE, include.missings = FALSE, custom.coef.names = c("Coup", "Ceasefire", "GDP (logged)", "Level of Democracy", "Intensity Level", "Conflict Duration (logged)", "Relative Strengths", "Territorial Control", "Multirebel Conflict", "Territorial Conflict", "Failed Coup", "Successful Coup"))

pa5 = coxph(Surv(istart_pa, ifinish_pa, pa) ~ military_coup + all_cf + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terr + strata(enum_pa) + frailty(gwno_a,distribution="gamma",method="em"), method=c("efron"), data = pa)

pa6 = coxph(Surv(istart_pa, ifinish_pa, pa) ~ military_coup + all_cf + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terr + multirebel + territorial + strata(enum_pa) + frailty(gwno_a,distribution="gamma",method="em"), method=c("efron"), data = pa)


pa7 = coxph(Surv(istart_pa, ifinish_pa, pa) ~ catcoupm + all_cf + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terr + strata(enum_pa) + frailty(gwno_a,distribution="gamma",method="em"), method=c("efron"), data = pa)

pa8 = coxph(Surv(istart_pa, ifinish_pa, pa) ~ catcoupm + all_cf + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terr + multirebel + territorial + strata(enum_pa) + frailty(gwno_a,distribution="gamma",method="em"), method=c("efron"), data = pa)

# hazard ratios
pa5h = exp(coef(pa5))
pa6h = exp(coef(pa6))
pa7h =  exp(coef(pa7))
pa8h = exp(coef(pa8))


texreg(list(pa5, pa6, pa7, pa8), 
       override.coef = list(pa5h, pa6h, pa7h, pa8h), include.aic = FALSE, include.rsquared = FALSE, include.zph = FALSE, include.missings = FALSE, custom.coef.names = c("Military Coup", "Ceasefire", "GDP (logged)", "Level of Democracy", "Intensity Level", "Conflict Duration (logged)", "Relative Strengths", "Territorial Control", "Multirebel Conflict", "Territorial Conflict", "Failed Military Coup", "Successful Military Coup"))
       


## summary statistics
sum_pa = pa %>% mutate(log_gdp = log(gdp), log_time = log(time)) %>%  select(pa, coup, military_coup, all_cf, log_gdp, v2x_polyarchy, intensity_level, log_time, rebstrength_num, terr, multirebel, territorial)

stargazer(as.data.frame(sum_pa), covariate.labels = c("Peace Agreement", "Coup", "Military Coup", "Ceasefire", "GDP (logged)", "Level of Democracy", "Intensity Level", "Conflict Duration (logged)", "Relative Strengths", "Territorial Control", "Multirebel Conflict", "Territorial Conflict"))
```




```{r}
m.pa1 = matchit(coup ~ all_cf + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terr + multirebel + territorial, data = pa, method = "full", caliper = 0.05, std.caliper = FALSE)
matched_pa1 = match.data(m.pa1)

pa_matching1 = coxph(Surv(istart_pa, ifinish_pa, pa) ~ coup + all_cf + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terr + multirebel + territorial, weights = weights, 
      cluster = gwno_a, data = matched_pa1)



pa_matching2 = coxph(Surv(istart_pa, ifinish_pa, pa) ~ coup + all_cf + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terr + multirebel + territorial, weights = weights, 
      cluster = subclass, data = matched_pa1)





m.pa2 = matchit(coup ~ all_cf + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terr + multirebel + territorial, data = pa, method = "full", caliper = 0.01, std.caliper = FALSE)
matched_pa2 = match.data(m.pa2)

pa_matching3 = coxph(Surv(istart_pa, ifinish_pa, pa) ~ coup + all_cf + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terr + multirebel + territorial, weights = weights, 
      cluster = gwno_a, data = matched_pa2)



pa_matching4 = coxph(Surv(istart_pa, ifinish_pa, pa) ~ coup + all_cf + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terr + multirebel + territorial, weights = weights, 
      cluster = subclass, data = matched_pa2)

new.names.pa = c(
  "coup" = "Coup", 
  "all_cf" = "Ceasefire",
  "log(gdp)" = "GDP (logged)", 
  "v2x_polyarchy" = "Level of Democracy",
  "intensity_level_2" = "Intensity Level",
  "log(time)" = "Conflict Duration (logged)", 
  "rebstrength_num" = "Relative Strengths",
  "terr" = "Territorial Control",
  "multirebel" = "Multirebel Conflict",
  "territorial" = "Territorial Conflict",
  "distance" = "Distance"
)




pa1h =  exp(coef(pa_matching1))
pa2h =  exp(coef(pa_matching2))

pa3h =  exp(coef(pa_matching3))
pa4h =  exp(coef(pa_matching4))



texreg(list(pa_matching1, pa_matching2, pa_matching3, pa_matching4), include.aic = FALSE, include.rsquared = FALSE, include.zph = FALSE, include.missings = FALSE, override.coef = list(pa1h, pa2h, pa3h, pa4h), custom.coef.names = c("Coup", "Ceasefires", "GDP (logged)", "Level of Democracy", "Intensity Level", "Conflict Duration (logged)", "Relative Strengths", "Territorial Control", "Multirebel Conflict", "Territorial Conflict"))

love.plot(m.pa1, var.names = new.names.pa)
love.plot(m.pa2, var.names = new.names.pa)



print(xtable(as.data.frame(bal.tab(m.pa1, var.names = new.names.pa)$Balance), digits = 3), sanitize.rownames.function = function(x) new.names.pa[x])

print(xtable(as.data.frame(bal.tab(m.pa2, var.names = new.names.pa)$Balance), digits = 3), sanitize.rownames.function = function(x) new.names.pa[x])


```




```{r}
### Visualization (Peace Agreement)
mod_pa = coxph(Surv(istart_pa, ifinish_pa, pa) ~ coup + all_cf + log(gdp) + v2x_polyarchy+  intensity_level +  log(time) +  rebstrength_num +terr + multirebel + territorial, data = pa)



new_pa = data.frame(
  coup = c(0, 1),  # Example: Compare coup = 0 vs. coup = 1
  all_cf = mean(pa$all_cf, na.rm = TRUE),
  gdp = mean(pa$gdp, na.rm = TRUE),
  v2x_polyarchy = mean(pa$v2x_polyarchy, na.rm = TRUE),
  intensity_level = mean(pa$intensity_level, na.rm = TRUE),
  time = mean(pa$time, na.rm = TRUE),
  terr = mean(pa$terr, na.rm = TRUE),
  rebstrength_num = mean(pa$rebstrength_num, na.rm = TRUE),
  multirebel = mean(pa$multirebel, na.rm = TRUE),
  territorial = mean(pa$territorial, na.rm = TRUE)
)


surv_pa = survfit(mod_pa, newdata = new_pa)

ggsurvplot(surv_pa, data = pa, 
           conf.int = FALSE,   
           legend.title = "", 
           legend.labs = c("No Coup", "Coup"),
           title = "",
           xlab = "Time (Years)",
           ylab = "Survival Probability",
           ggtheme = theme_minimal(),
           xlim = c(0,20),
           palette = c("lightgrey", "darkgray"),  
           linetype = c("dotted", "solid"))
```




```{r}
cf_yearly = cf %>%
  group_by(year) %>%
  summarise(ceasefire_count = sum(all_cf == 1, na.rm = TRUE)) %>%
  ungroup()


pa_yearly = pa %>%
  group_by(year) %>%
  summarise(peace_agreement_count = sum(pa == 1, na.rm = TRUE)) %>%
  ungroup()

cf_pa_yearly = full_join(cf_yearly, pa_yearly, by = "year") %>%
  replace_na(list(ceasefire_count = 0, peace_agreement_count = 0))


cf_pa_yearly_long = cf_pa_yearly %>%
  pivot_longer(cols = c(ceasefire_count, peace_agreement_count),
               names_to = "event_type", values_to = "count") %>%
   mutate(event_type = case_when(
    event_type == "ceasefire_count" ~ "Ceasefires",
    event_type == "peace_agreement_count" ~ "Peace Agreements"
  ))

ggplot(cf_pa_yearly_long, aes(x = year, y = count, color = event_type)) +
  geom_line(size = 0.7) +
  geom_point() +
  labs(title = "",
       x = "Year", y = "Count",
       color = "") +
  theme_minimal()
```
